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Abstract 

Maxwell's equations are cast in the form of the Schrodinger equation. The Lanc- 
zos propagation method is used in combination with the fast Fourier pseudospectral 
method to solve the initial value problem. As a result, a time-domain, unconditionally 
stable, and highly efficient numerical algorithm is obtained for the propagation and 
scattering of broad-band electromagnetic pulses in dispersive and absorbing media. As 
compared to conventional finite-difference time-domain methods, an important advan- 
tage of the proposed algorithm is a dynamical control of accuracy: Variable time steps 
or variable computational costs per time step with error control are possible. The 
method is illustrated with numerical simulations of extraordinary transmission and 
reflection in metal and dielectric gratings with rectangular and cylindrical geometry. 
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I Introduction 



Recent developments in photonics and nanostructure materials [1, 2] have increased inter- 
est in efficient and accurate algorithms for numerical simulations of the propagation and 
scattering of short (broad band) laser pulses in generic passive (dispersive and absorbing) 
media. Time-domain approaches for solving the Maxwell's equations might be more suit- 
able for this purpose than frequency domain methods because the sought-for information, 
e.g., the scattering matrix, can be obtained within a desired frequency range by a single 
propagation. Coupled with laser ellipsometry of broad band pulses, fast simulations of ex- 
pected resonance patterns in the scattering amplitude appear to be an efficient tool to control 
quality of manufactured photonic devices. Unconditionally stable algorithms are especially 
advantageous for such tasks because of their applicability to practically all materials and 
geometries without any assessment of admissible values of the system parameters. Another 
attractive property of time domain methods is their universality. The very same algorithms 
can be used to calculate static properties of the system (e.g., a band structure of photonic 
crystals), to simulate the electromagnetic pulse propagation in non-linear materials as well 
as in media with time-dependent properties. 

The advantages of time-domain methods have been for a long time recognized in quantum 
mechanics where they are extensively used in the fields of chemical reaction dynamics [3], laser 
- matter interactions [4] , etc. Highly efficient and accurate tools have been developed for the 
wave packet propagation and analysis of the results [3, 5, 6, 7, 8, 9]. Since Maxwell's equations 
can be cast in the form of the Schrodinger equation, it is then natural to extend time- 
domain methods of quantum mechanics to numerical electrodynamics. Some realizations 
of this idea are rooted to the path integral representation of quantum theory (the Lie- 
Trotter product formula [10] or the split operator method [5, 11, 12]). The others exploit 
polynomial approximations of the fundamental solution of the Schodinger equation. For 
instance, the Chebychev time-propagation technique has been recently used to simulate the 
electromagnetic pulse propagation in non- absorbing media [13, 14]. 

Here it is proposed to use the Lanczos algorithm [15] to obtain an unconditionally stable, 
time-domain solver of Maxwell's equations for passive media. The method allows for a 
dynamical control of accuracy, meaning that computational costs are constantly optimized 
in due course of simulations with error control. In brief, the approach can be summarized 
as follows. Maxwell's equations are written in the form of the Schrodinger equation which is 
then solved by the Lanczos propagation scheme [16, 5] (Section II). The difference with the 
well studied quantum mechanical case is that the wave function is a multi-dimensional vector 
field and the Hamiltonian is non-Hermitian for absorbing media. The split operator method 
[5, 11, 12] has been used to include attenuation into the Lanczos propagation scheme, while 
preserving its unconditional stability (Section III). The action of the Hamiltonian on the 
wave function is computed by means of the Fourier pseudospectral method introduced in 
[17]. 

The accuracy of the method is investigated in Section IV. In Section V the Lanczos 
propagation scheme is applied to various gratings. In particular, a resonant extraordinary 
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reflection of a periodic array of parallel dielectric cylinders is observed. This effect is similar to 
the Wood anomalies [18] and related to the existence of stationary (trapped) electromagnetic 
waves with wave vectors parallel to the discrete translation symmetry axis of the system. 
Simulations of the scattering of broad band pulses on metallic grating and grooves, whose 
dielectric properties are described by the Drude model, are performed to demonstrate that 
the Lanczos propagation scheme is able to reproduce the results known in the literature 
and obtained by different means (by finite differencing schemes or by the scattering matrix 
method). 



II The Lanczos method for Maxwell's equations 

Consider first the case of non-dispersive media. Let D and B be electric and magnetic 
inductions, respectively, and E and H the corresponding fields so that D = eE and B = /xH 
where e and fi are positive, symmetric, position dependent matrices for generic non-isotropic 
and non-homogeneous media. For isotropic media, e and fi are scalars. At interfaces of 
different media, the boundary conditions are enforced dynamically, that is, e and /i are 
allowed to have discontinuities. Maxwell's equations are rewritten as: 

The over-dot denotes the time derivative, and c is the speed of light in the vacuum. One can 
also use the electromagnetic inductions as independent variables instead of the fields. The 
Hamiltonian TL must then be modified accordingly. The initial value problem is solved by 
applying the evolution operator (or the fundamental solution) to the initial wave function 

^(t) = e _<w V(0) • (2.2) 

The Hamiltonian is a Hermitian operator, TO = TL, with respect to the measure scalar 
product 

(V>i, $2) = J (Di ' E 2 + B 1 ■ H 2 ) dr = J Vl/# 2 dr . (2.3) 

The symmetric positive matrix k is block-diagonal, with the blocks being e and fi. The 
norm of the wave function with respect to the scalar product (2.3) is proportional to the 
electromagnetic energy and is conserved because the evolution operator is unitary. 

In numerical simulations, the Hilbert space is projected onto a finite dimensional Eu- 
clidean space so that tp becomes a vector whose components are values of the wave function 
at sites of a finite spatial grid. In the grid representation, TL is a matrix. If the Hamiltonian 
is Hermitian, it is then convenient to have TL as an explicitly Hermitian matrix. In the 
Maxwell theory, this can be achieved if, before projecting onto the grid, the wave function 
and the Hamiltonian are scaled 

$ _> K ~V^ ^ n K -V2 W/c i/2 ( 2 _4) 
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In the representation (2.1) we have 



/ E , H, ^_ ityx -i/2 Vxe -i/ 2 o J • (2-5) 

The scaled Hamiltonian is Hermitian with respect to the conventional scalar product in the 
space of square integrable functions, and, hence, it is a Hermitian matrix, when projected 
onto the grid. The action of spatial derivatives is calculated within the pseudospectral 
approach based on the Fourier grid representation of the wavefunction and the fast Fourier 
transform. In what follows, only consecutive actions of 7i on wave functions are required. 

A direct use of (2.2) implies a diagonalization of H, which is not feasible if the matrix 
size is too large. Various numerical approximations are based on the semigroup property of 
the evolution operator 

^(t + At) = e- iAtn ^(t) . (2.6) 

In a local propagation scheme the exponential can be approximated by a polynomial for a 
sufficiently small time step At. The basic idea of the Lanczos propagation method is that the 
exact solution ip(t + At) is projected onto the Krylov subspace associated with the initial 
state ip(t) and the Hamiltonian, ip(t + At) -> ip^(t + At) = V n ip(t + At) G K n , where 
V\= V n , V 2 n = V n , and 

K n = Span (^(t),^(t),..,H n -VW) • 
The accuracy of such an approximation is 0(At n ). The Hamiltonian is projected accordingly, 

U -> «W = V n HV n . Thus, 

i){t + At) w ^ n \t + At) = e- <Atw( "V (n) (*) • ( 2 - 7 ) 

The projection is done via an orthonormal basis for K n which is constructed by means 
of the Lanczos recursion algorithm [15]. In this basis, the matrix is Hermitian and 
tridiagonal. Typically, just a few orders are sufficient (n < 9) so that n is much smaller than 
the dimension of TC and the matrix 7iW can easily be diagonalized. The dimension n may 
be set differently at each time step, depending on the current vector ip(t), and is determined 
by a pre-set required accuracy. In particular, it allows to avoid excessive actions of Tt on the 
wave function. This feature leads to a dynamical optimization of computational costs with 
error control, which is one the greatest advantages of the Lanczos method. 

A detailed discussion of the Lanczos recursion algorithm and its application to the wave 
packet propagation can be found elsewhere [16, 5]. Here only a brief summary is given with 
notations used later in the text. Let ip = ip{t) where ip{t) is assumed to be normalized so 
that H^oll — 1- Due to the linearity of the Schrodinger equation one can always scale ip by 
a number and rescale it back after applying the infinitesimal evolution operator. Define 

a = (^ ,^o) , (2.8) 
0x = (W-aoM), (2-9) 

V>i = h/WM ■ (2.10) 
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By construction, ipi and ipo are orthonormal. For k — 2, 3, n — 1 the rest of the basis for 



K n is generated by the recursion relation 

a fc _i = (V>fc-i,7#fc-i) , (2-H) 

/3 fc _ 2 = (ip k - 2 ,Hipk-i) , (2.12) 

fc = (H - a fc _i)^ fc _i - p k ^k-2 , (2.13) 

= fc /||0 fe || • (2.14) 



By construction, the vector is a linear combination of V'j-ij V'jj an d V'j+i- Hence, in 
the Lanczos basis the matrix = (ipi^Tiipj) is tridiagonal. Elementary calculations show 
that the diagonal elements are Tiff = Oij = aj, the upper and lower superdiagonals are 

Let C/ be a unitary transformation such that WH^U is a diagonal matrix, and Ej be 
eigenvalues of HS n \ The approximate solution (2.7) is obtained by expanding the wave 
function over the Lanczos basis and solving the Schodinger equation for the expansion coef- 
ficients: 

n— 1 n—1 

^ n \t + At) = Ujk e- iAtE * U j0 ^ = J2 c ^t)^k , (2.15) 

k,j=0 fc=0 

where the initial condition c&(0) = 5ko has been taken into account. Since TiS 7 ^ is Hermitian, 
the evolution preserves the norm 

||^W(t + A0|| 2 = ||^ B) (0ir = IW 2 = l- (2-16) 
Thus, the algorithm is unconditionally stable because the norm of the amplification matrix 
G in) (At), defined by ^ n \t + At) = g {n \At)^{t), is uniformly bounded, ||£ (n) (£)|| < 1, for 
all parameters of the Hamiltonian and At > 0. 

The accuracy of the algorithm can be estimated from the following observation [16]. The 
norm of a projection of the exact solution onto the orthogonal complement of K n can be 
used as a measure of accuracy of the Lanczos algorithm. By expanding the exponential in 
the right hand side of (2.6) into the Taylor series, it is clear that the contribution of the term 
(AtT-C) n+1 ip(t), which has no projection onto K n , can only be captured by the approximate 
solution if the larger Krylov space JC n+ 2 is used in the Lanczos algorithm, which, in turn, 
implies that the vector Cj(At) acquires two additional components. Thus, the accuracy of 
the Lanczos algorithm can be controlled, for example, by demanding that the absolute value 
of c n _i(At) is less than a specified small number e. Note that |c n _i(At)| ~ 0{At n ~ 2 ) as one 
can deduce from (2.15) and the tridiagonal structure of in the Lanzcos basis. To ensure 
that the norm of the projection of ip(t + At) onto the orthogonal complement of K n is small, 
we demand that 

|c„_ 3 (At)| 2 + |c„_ 2 (At)| 2 + |c n _!(At)| 2 < e (2.17) 

where e ~ 1CT 14 in our calculations. To satisfy (2.17), the time step At, or the dimension of 
the Krylov subspace n, or both can be varied to minimize computational costs. This is the 
aforementioned dynamical control of accuracy in the Lanczos propagation method. In our 
simulations, At has been kept fixed, while (2.17) has been used to determine a minimal n 
for each time step. 
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Ill Including attenuation by the split method 



The response function of a passive medium in an applied electromagnetic field must satisfy 
the causality condition. A common way to model the causal response function is to assume 
that the medium polarization and magnetization satisfy a linear differential equation in time 
in which a non-homogeneous term is proportional to the applied field (for linear media). The 
Maxwell's equations in passive media appear then to be a system of (high-order) differential 
equations to which numerical algorithms are applied [19, 20]. Any system of high-order 
differential equations can be converted into a system of first-order differential equations by 
introducing auxiliary dynamical variables. This idea is used to convert Maxwell's equations 
for passive media into the Schrodinger equation (2.1) in which the wave function contains ad- 
ditional components that describe dynamics of the medium polarization and magnetization. 
Due to absorption the time evolution is no longer unitary. 

It must be noted that absorption of the wave packet is required in numerical simulations 
of scattering problems in which the pulse shape is to be computed in the asymptotic region. 
Indeed, when the front edge of the pulse reaches the grid boundary, it will be reflected or 
re-appear on the other side of the grid, depending on the boundary conditions. To avoid 
an artificial interference of the scattered pulse with itself, a layer of an absorbing medium is 
necessary at the grid boundary [21]. 

Here a simple way is proposed to include the attenuation of the wave packet amplitude 
into the Lanczos method, while maintaining the unconditional stability of the algorithm. 
The procedure is illustrated with the Drude model of metals. 

Let Ti = Ti,o — iV where Ti^ = TCo and = V. The system is absorbing and, therefore, 
V must be a positive semidefinite operator, that is, for any ip, (■?/>, Vip) > 0. This readily 
follows from the condition that the norm of a solution of (2.1) cannot increase with time. 
The exact time evolution (2.6) is approximated by means of the Lie- Trotter formula 

ij(t + At) = e~ Atv / 2 e- lAtHo e~ Atv / 2 ^(t) + 0(At 3 ) . (3.1) 

The action of the exponential of 7Yo is computed by the Lanczos method as before. The 
attenuation potential V typically does not involve spatial derivatives and, hence, the action 
of its exponential on a wave function is far less expensive than that for Ho. The norm 
of any power of the amplification matrix still remains uniformly bounded by one because 
|| e -Atv/2|| <- i £ or ^ > q H ence th e unconditional stability is preserved. 

Let us turn to the Drude model which is used in numerical simulations presented below. 
Another popular model, a multi-resonance Lorenz model, can be treated similarly. Let 
D = E + P and B = H. In the Drude model, the medium polarization is described by the 
second order differential equation 

P + ??P = ^E, (3.2) 

where rj > is the attenuation constant and u p is the plasma frequency. Equation (3.2) 
must be solved with zero initial conditions, P(0) = P(0) = 0. Define an auxiliary field Q 
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by P = uj p Q. Rewriting the Maxwell's equations and (3.2) in terms E, B, and Q and their 
first-order time derivatives, the Schrodinger equation is obtained in which 

(E\ / icVx -iu p \ 

B , ft = I -icVx . (3.3) 
Q / \ iuj p -ir\ j 

The Hamiltonian is Hermitian when 77 = (no attenuation). The attenuation potential V is 
a diagonal matrix, diag (0, 0, 77), that is positively semidefinite since 77 > 0. 

For an absorber at the grid boundaries, a layer of a conducting medium has been used 
with a position dependent conductivity a. As the induced current in a conducting medium 
has the form J = crE, the matrix V is changed to diag (— Ana, 0, 77). The function a is 
constructed according to the frequency band of the initial pulse. 



IV Free space propagation. Phase and amplitude er- 
rors 

To illustrate the efficiency of the Lanczos time-propagation method, we compare it with a 
widely adopted Second Order Finite Differencing (SOD, or leapfrog) propagation method 
[5, 20, 22], using the simplest example of the electromagnetic pulse propagation in vacuum. 
The action of the Hamiltonian on wave functions in the Lanczos and leapfrog methods are 
done in the same way, that is, by the fast Fourier pseudospectral method on the same grid. 

Consider a Gaussian wave packet linearly polarized along the y axis and propagating 
along the z axis. The amplitude of the fields at the initial time t — is given by 

E y {z) = e -* 2 l D2 e ik ° z , H x (z) = -E y {z) , (4.1) 

where k = 5.5/ D, and D determines the width of the wave packet. The carrier wave length 
A = 27r//c so that D = 0.875A. We take D = 1.75 [im, or A = 2 fiin. The step of the grid 
is Az = 0.1D. An exact solution directly follows from (4.1) E y (z,t) = E y (z — ct). The 
wave packet propagates in the direction of positive z. With our settings the pulse duration 
is about 25 fs. 

Numerical solutions are obtained by the Lanczos and leapfrog algorithms for the Schro- 
dinger equation (2.1) in which e — fx — 1. Recall that the leapfrog propagation scheme is 
based on the third-order finite difference approximation of the time derivative 

ip(t + At) = ip(t - At) - 2iAtHip{t) . (4.2) 

The scheme is conditionally stable, and the time step must be chosen accordingly. The 
simulated electric field is recorded by a detector placed at z — z^et = 18D. Its phase and 
amplitude are compared with those of the exact solution. For a signal E(t) = E (t)e l(p ^\ 
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where E (t) = \E(t)\, the phase and amplitude errors are defined, respectively, by 

Yp^xact ^approx | ^£jexact j^approx | 

^ ~ (pexact ' ^ Jjjexact ' (^-3) 

The errors 5 P ' A are plotted respectively in Figs 1 and 2 as functions of S — (z dct — ci) / D, 
the position of the pulse center relative to the detector measured in units of D. The results 
are shown for \S\ < 2.5 where the signal on the detector is sufficient. Dashed and solid lines 
correspond to the leapfrog and Lanczos methods, respectively, for various settings of the 
time step. 

The time step for the black dashed line is a reference time step, At xs 0.01 fs. If N H is 
the number of elementary operations required to compute the action of the Hamiltonian on a 
wave function, then the total number of operations reads N = sNnN t , where s is the number 
of actions of the Hamiltonian per a time step, N t = t/ At is the total number of time steps. 
For the leapfrog method, s = 1 for all time steps. In the Lanczos method, s — n — 1, with n 
being the dimension of the Krylov space. Despite that the dynamic control of accuracy has 
been activated, de facto n does not vary in due course of simulations in vacuum. 

Let N — N for the black dashed curve. The red dashed curve is obtained by reducing 
the time step, At = At /2, and, hence, the total number of operations increases accordingly, 
iV = 2Nq. In the Lanczos method, the black solid curve corresponds to At = lOAto and 
s = 7, the blue solid curve to At = 5At and s = 7, and the red one to At = 2.5At 
and s = 6. The total number of operations is, respectively, iV = 0.7A^ , iV = 1.4A^ , and 
iV = 2.4iV . It is readily seen that at roughly the same number of operations, the Lanczos 
algorithm has phase and amplitude errors that are less than those in the leapfrog method 
by several orders of magnitude. 

A few remarks are in order. There are, of course, algorithms that would be more effi- 
cient than the Lanczos propagation method in free space. For instance, the split propagation 
method [12] essentially reproduces an exact solution and is also unconditionally stable. How- 
ever, the split method would not be applicable when the Hamiltonian involves products of 
operators that depend on spatial derivatives and positions. The accuracy of the leapfrog 
scheme can be improved by, for example, taking into account the next term of the Taylor 
expansion oiip(t± At) in powers of At in (4.2) [23], 

-2iAtHip -> -2iAtH(l - At 2 H 2 /3)ip . 

In this case, s = 3. The method is still conditionally stable where the stability condition of 
the SOD, At\\n\\ < 1, changes accordingly to At\\H(l - At 2 H 2 /3)\\ < 1. Even though s 
has tripled, the new stability condition allows one to increase the time step by the factor of 
2.1. Therefore the total number of operations increases only slightly. The accuracy of the 
scheme will be of 0(At 5 ) which is still not as high as in the Lanczos method, 0(At n ) with 
n = 8, 7 in the above examples. 

Some care should be taken regarding a known drawback of the Lanczos algorithm - a 
possible loss of orthogonality of basis functions due to round-off errors [15, 24]. This is why 
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the time step has to be adjusted so that only low dimensional Krylov spaces, n < 9, are 
invoked in contrast to the conventional use of the Lanczos method for solving linear systems. 



V Applications to gratings 

In this section the Lanczos propagation scheme is applied to the scattering of a broad band 
wave packet on nanostructure periodic materials such as gratings and grooves. We are 
particularly interested in transmission (reflection) properties currently being a subject of 
intense research [2, 26, 27, 28]. The results obtained here are compared with those available 
in the literature. The time-dependent approach allows us to underline the role played by 
trapped modes or resonances in the existence of extraordinary transmittance and reflectance 
of periodic structures. The longer lives a trapped mode, the more narrow resonance occurs 
in the reflection and/or transmission coefficient. 

All systems considered here have a translation symmetry along one of the Euclidean axes, 
chosen to be the y axis. The structures are periodic along the x axis with the period D g , 
and the z direction is transverse to the structure. The initial wave packet is Gaussian and 
propagates along the z axis. Its spectrum is broad enough to cover the frequency range of 
interest. The zero diffraction mode is studied for wavelengths A > D g so that reflected and 
transmitted beams propagate along the z-axis. As in our previous work [25] we use a change 
of variables to enhance the sampling efficiency in the vicinity of medium interfaces so that 
the boundary conditions at sharp interfaces are accurately reproduced by the Fourier-grid 
pseudospectral method. A typical size of the mesh corresponds to —15D g < z < 15D g , and 
— 0.5D g < x < 0.5D g with, respectively, 512 and 128 knots. The frequency resolved trans- 
mission and reflection coefficients are obtained via the time-to-frequency Fourier transform 
of the signal on "virtual detectors" placed at some distance in front and behind the slab 
with a periodic structure [29]. 

V.l Array of dielectric cylinders 

The significance of trapped modes is first illustrated with a periodic array of non- dispersive 
dielectric cylinders, the system which has not received as much attention as metal or dielectric 
gratings. Consider an array of parallel, periodically positioned, dielectric cylinders in vacuum 
oriented along the y axis. The radius R of cylinders is small as compared to the array period 
D g = 1.75/j.m. In simulations, the ratio R/D g is taken to be 0.0857. The incident wave 
packet is linearly polarized. The electric field is oriented along the y axis, i.e., parallel to the 
cylinders (the so called TM polarization). The Hamiltonian for the Lanczos scheme has the 
form (2.1) where fi = 1. 

In Fig. 3 the reflection coefficient 1Z is shown as a function of the wave length expressed 
in units of D g . In the Schrodinger formulation of Maxwell's theory the norm of the wave 
function is proportional to the total electromagnetic energy. Hence, for a lossless medium the 
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transmission T can simply be obtained from the energy conservation: T +1Z = 1. Recall that 
the Lanczos propagation method preserves the norm. The solid-blue and dashed-red curves 
correspond, respectively, to e = 2 and e — 4. As one can see the array becomes a perfect 
reflector within a fairly narrow wavelength range centered at the resonant wavelength that is 
slightly larger than the period D g . Similar results have been obtained for dielectric grating 
structures. The resonant pattern is associated with the so-called Wood anomalies [18], and 
can be explained by the existence of trapped modes or guided wave resonances [25, 30]. The 
widths of the resonances in the reflection (transmission) coefficient are determined by the 
lifetime of a corresponding quasi-stationary trapped mode which is a standing wave along 
the x axis and is excited by the incoming wave. 

The existence of trapped modes can easily be inferred from the temporal evolution of 
the electromagnetic field. Figure 4 shows the transmitted electric field as a function of time 
measured by a detector placed behind the layer of dielectric cylinders. The main transmitted 
pulse is clearly visible. It has a significant amplitude and duration about 25 fs. After the 
main pulse passes the array, it leaves behind an excited quasi-stationary mode which looses its 
energy by radiating almost monochromatic waves with the same amplitude, but an opposite 
phase, in the transmission and reflection directions. The lasing effect of the trapped mode 
appears as exponentially dumped oscillations coming after the main signal. The exponential 
decay due to a finite lifetime of the quasi-stationary state is clearly seen. By the symmetry, 
the same lasing effect is registered by a detector placed in front of the layer (not shown here). 
A 100% reflection at the resonant frequency can be understood from the fact that the field 
emitted by the trapped mode in the transmission direction and the corresponding frequency 
component of the initially transmitted pulse have an opposite phase, thus compensating each 
other. The solid-blue and dashed-red curves correspond, respectively, to e = 2 and e — 4. 
The radiation coming from the narrow resonance (the blue curve) has a lower amplitude and 
a much longer duration. The lifetime of the trapped mode in this case is in the picosecond 
range, i.e., thousand times longer than the initial pulse duration. Note that the more narrow 
resonance is the less energy gets trapped from the initial pulse. This explains the amplitude 
difference of the blue and red curves. Finally, the concept of trapped modes localized on 
successive layers and interacting with each other provides a theoretical framework for the 
light propagation in layered structures such as photonic crystal slabs [31]. 



V.2 Metal gratings and grooves 

Metal gratings and grooves have been extensively studied in micro- wave and optical domains 
[26, 27]. The purpose of this section is to show that the Lanczos propagation method can 
successfully be applied to metals described by the Drude model. The Hamiltonian has the 
form (3.3). The attenuation and the plasma frequency are taken to be representative for 
silver: uo p = 9eV, and rj = O.leV [26]. The grating geometry is sketched in the inset of Fig. 
5. The grating period is D g = 1.75 //to, the thickness (along the z axis) is h — 0.8 //to, and 
the grating width a = 0.3 //to. The corresponding grooves are obtained by attaching a solid 
metal plate on one side of the gratings so that no transmission is possible. The polarization 
of the incident wave packet is such that the electric field vector is oriented along the x axis, 
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i.e., perpendicular to the gratings (the so called TE polarization). The difference with the 
non-dispersive case discussed above is the presence of attenuation. The trapped mode looses 
its energy due to (non-perfect) conductivity of the metal. This leads to broadening of the 
resonance. 

In Fig. 5 the dashed red and solid blue curves represent the transmission and reflection 
coefficients, respectively, as functions of the wavelength expressed in units of the grating 
period, D g . The resonance is again associated with the existence of a trapped stationary 
wave in the grating. The transmittance does not reach 100% due to dissipative loss of energy 
in the Drude metal. While for a lossless medium the sum of the reflection and transmission 
coefficients must be one, this is not the case for the Drude metal (the dashed-dotted green 
curve in Fig. 5). The maximal loss of energy corresponds to the resonant wavelength. It is 
easily understood because the trapped mode remains in contact with the metal much longer 
than the main pulse, and, therefore, can dissipate more energy through exciting surface 
electrical currents in metal. The black curve in Fig. 5 shows the reflectance of the grooves. 
Since the light cannot be transmitted through the grooves, a resonance structure in the 
reflection coefficient is directly related to the enhanced energy loss at the wavelength of the 
trapped mode. Note that as compared to the metal gratings, the resonance is broadened 
and shifted to the lower frequencies (larger wavelength). The results obtained here are in a 
full agreement with previous theoretical and numerical analysis [26, 27]. 



VI Conclusions 



It has been demonstrated that the Lanczos algorithm can be used to develop a highly efficient, 
accurate, and unconditionally stable propagation scheme to simulate scattering of broad band 
electromagnetic pulses in passive media. The accuracy and efficiency of the algorithm have 
been illustrated with an example of the electromagnetic wave propagation in vacuum. At 
the same computational costs, a significant reduction of phase and amplitude errors has 
been observed in the Lanczos propagation method as compared to the second-order finite- 
difference (leapfrog) scheme. 

As an example of possible applications, the Lanczos propagation method has been applied 
to study resonant transmission and reflection of various periodic nano-structures: An array 
of periodically placed parallel cylinders made of a non-dispersive dielectric material, metal- 
lic gratings and grooves. The time-domain study clearly demonstrates the role played by 
quasi-stationary (trapped) electromagnetic waves supported by the corresponding periodic 
structure in the extraordinary transmission (reflection) properties of the grating. The re- 
sults for metallic gratings and grooves coincide with those obtained earlier by means of other 
numerical algorithms and are also in agreement with theoretical studies. The unconditional 
stability of the Lanczos propagation scheme for media with attenuation has been achieved 
via the split method, which reduces the accuracy. It is possible to restore the accuracy up 
to the level gained for non-absorbing media. However, stability conditions require a further 
study that will be reported elsewhere. 
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In summary, the Lanczos algorithm has been shown to lead to a highly accurate, efficient, 
and unconditionally stable time-propagation numerical solver for the Maxwell's equations. 
Variable time steps and/or variable computational costs with accuracy control are possible. 
The method is applicable to various electromagnetic systems (no restrictions on the Hamil- 
tonian). All these virtues are hardly available in other unconditionally stable algorithms in 
numerical electrodynamics of passive media. 
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Figure captions 

Fig. I. The phase errors for the propagation of an electromagnetic gaussian pulse 
in vacuum. Results are presented as a function of the position of the pulse center relative 
to the detector, S, measured in units of the pulse width D. Dashed and solid curves cor- 
respond, respectively, to the leapfrog and Lanczos propagation methods. Different colors 
represent computational costs of simulations measured as the total number of actions of the 
Hamiltonian on the wave function for fixed propagation time. Further details are given in 
the text. 

Fig. 2. The amplitude errors for the propagation of an electromagnetic gaussian 
pulse in vacuum. Results are presented as a function of the position of the pulse center 
relative to the detector, S, measured in units of the pulse width D. Dashed and solid curves 
correspond, respectively, to the leapfrog and Lanczos propagation methods. Different colors 
represent computational costs of simulations measured as the total number of actions of the 
Hamiltonian on the wave function for fixed propagation time. Further details are given in 
the text. 

Fig. 3 Calculated zero-order reflection coefficient for a periodic array of dielectric 
cylinders in vacuum described in the text. Results are presented as a function of the wave- 
length of the incident radiation measured in units of the period D g . The solid blue and dashed 
red curves correspond, respectively, to the array of cylinders with dielectric constants e = 2 
and e — 4. 

Fig. 4 The electric field measured by a detector placed behind the periodic layer 
of dielectric cylinders. Only the field corresponding to the zero-order transmitted wave 
propagating along the z-axis is represented. It is obtained by the Fourier analysis of the 
^-coordinate dependence of the field at the detector position. The signal is shown as a 
function of time measured in femtoseconds. The solid blue and dashed red curves correspond, 
respectively, to the array of cylinders with dielectric constants e = 2 and e — 4. 

Fig. 5. Calculated zero-order reflection and transmission coefficients for metallic grat- 
ings and groves described in the text. Results are presented as a function of the wavelength 
of the incident radiation measured in units of the period D g . The inset of the figure gives a 
schematic view on the grating geometry. The black line shows the reflection coefficient for 
metallic grooves. Blue and (dashed red) line shows the reflection (transmission) coefficient 
for metallic gratings. The sum of the reflection and transmission coefficients for metallic 
gratings is shown as the dashed-dotted green curve. Its deviation from 1 represents the loss 
of electromagnetic energy because of the absorption in metal. 
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